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Abstract 

We simulated the 3-band Hubbard model using the Constrained Path Monte 
Carlo (CPMC) method in search for a possible superconducting ground state. 
The CPMC is a ground state method which is free of the exponential scaling 
of computing time with system size. We calculated the binding energy of a 
pair of holes for systems up to 6 x 4 unit cells. We also studied the pairing 
correlation functions versus distance for both the d-wave and extended s-wave 
channels in systems up to 6 x 6. We found that holes bind for a wide range 
of parameters and that the binding increased as the system size is increased. 
However, the pairing correlation functions decay quickly with distance. For 
the extended s channel, we found that as the Coulomb interaction Ud on 
the Cu sites is increased, the long-range part of the correlation functions is 
suppressed and fluctuates around zero. For the d^2_^y2 channel, we found that 
the correlations decay rapidly with distance towards a small positive value. 
However, this value becomes smaller as the interaction or the system size 
is increased. 
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PACS Numbers: 74.10.-z, 74.20.+V, 71.10.Fd, 71.10.-w, 02.70.-c 



I. INTRODUCTION 

The discovery of high temperature superconductivity in copper oxide materials generated 
considerable effort to find simple electronic models that exhibit a superconducting phase. 
Since the mechanism for superconductivity is believed to lie within the Cu02 planes, a variety 
of two-dimensional models have been proposed, the three principal ones being the (one- 
band) Hubbard, three-band Hubbard, and t-J models. Unfortunately, conclusive evidence 
for a superconducting phase has not been found in any of these models.0 

In this paper, we will report the results of a quantum Monte Carlo (QMC) study of 
the three-band Hubbard model. Of the three models, it has been less intensively studied, 
partially because of the general belief that its low energy excitation spectrum is similar to the 
other two. Indeed in the strong coupling limit, both the one-band and three-band Hubbard 
models have the t-J model as an approximate limit. However, there still remains some 
controversy on whether one-band models, like the Hubbard and t-J models, are adequate 
to describe the low energy physical properties of the cuprate superconductors. One of our 
objectives is to study the possible existence of superconductivity in the three-band model 
in regions where model parameters are physical as opposed to regions where asymptotic 
models are clearly more appropriate. 

The most solid information about possible superconductivity in the one-band Hubbard 
model has come from a series of QMC calculations. For instance, using a finite temperature 
QMC method. White et al.i studied possible superconductivity in the Hubbard model. They 
found an attractive effective pairing interaction in the d^-i-y-i and extended s-wave channels. 
Moreo and Scalapinoli subsequently found pairing correlations but also found that they did 
not increase when the lattice size was increased. Their result, suggesting the absence of long 
range order, was consistent with an earlier QMC study by Imada and Hatsugai.@ 

The results of almost all QMC calculations of the Hubbard model however, are limited 
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by the fermion sign problem to high temperatures and small system sizes. These limitations 
had left open the possibility that superconductivity still lurked at larger systems and lower 
temperatures. Recently, a new QMC method, the Constrained Path Monte Carlo (CPMC) 
method,i was developed to get around the sign problem. Using this method, Zhang et al.i 
calculated ground-state pairing correlation functions for the Hubbard model as a function 
of the distance and found that as the system size or the interaction strength is increased the 
magnitude of the long-range part of the pairing correlation functions vanished for both the 
and extended s-wave channels. Although this method produces a variational solution, 
its results, together with the null results from previous QMC studies, are very discouraging 
for finding superconductivity in the one-band Hubbard model. 

Past numerical work on the three-band model, has been less definitive about the ex- 
istence of superconductivity. The work has mainly consisted of exact diagonalization 
computations@~i, where calculations of hole binding energies was emphasized, and QMC 
computations0~0 where magnetic and pairing calculations were emphasized. The exact 
diagonalization studies unequivocally established that holes can bind. In a number of cases 
binding only occurred for unphysical parameter ranges, and because of small systems sizes, 
the existence and symmetry of the pairing was difficult to access. The QMC studies, prob- 
ably even more limited by the sign problem than those for the one-band model, established 
the existence of an extended s-wave and dj.2_y2 wave attractive pairing interaction. However, 
none of these QMC studies reported the behavior of the pairing correlation functions at large 
distances, which is the relevant quantity to establish the existence of long range order. 

In the work reported here, we removed the limitations caused by the sign problem by 
applying the CPMC method to the three-band Hubbard model. We computed the binding 
energy for holes near half-filling for systems up to 6 x 4 unit cells. We found a wide range 
of physically relevant parameters for which the holes bind. As will be reported in more 
detail, we found that larger systems in fact make hole binding easier. We also find that 
the inclusion of the nearest neighbors Coulomb repulsion Vpd is not necessary for binding 
as was concluded in several of the exact diagonalization studies. For similar systems sizes 



as the previous simulations, but at zero temperature, we looked at the pairing correlation 
functions as a function of distance and found that the pairing correlation functions decay 
quickly with distance and are suppressed with increasing lattice size and Ud- The extended 
s-wave channel is much more strongly suppressed than the d^^-y^ channel. These trends are 
similar to those found in QMC simulations of the one-band model. 

The remainder of our report is organized as follows: in Section |T| we define the Hamil- 
tonian and the physical quantities calculated and discuss the choice of model parameters. 



In Section |T| we briefly describe the CPMC method, and then in Section |^ we present the 
results for both the binding energy for holes and the pairing correlation functions. Finally 
in Section 0, we discuss in detail our main conclusions. 



II. THREE-BAND HUBBARD HAMILTONIAN 

As proposed by Emeryjil the three-band model physically mimics the CUO2 layer by 
having one Cu and two O atoms per unit cell, with the Cu atoms arranged on a square 
lattice and the O atoms centered on the edges of the square unit cells. In this layer, Emery 
assumed that the relevant orbitals are just those of copper 3d^2_y2 and oxygen 2px and 2py. 

The Hamiltonian has the form: 

<j,k>a ja j 

+edJ2nt + UdY.« (1) 

icr i 

+ Vpd ^i^^j + Y tpdi^tPja + P]adia) 

<i,j> <i,3>" 

In writing the Hamiltonian, we adopted the convention that the operator d\„ creates a hole 
at a Cu 2>dx2-y2 orbital and p]^ creates a hole in an O 2px or 2py orbital. Ud and Up are 
the Coulomb repulsions at the Cu and O sites respectively, and are the corresponding 
orbital energies, and Vpd is the nearest neighbor Coulomb repulsion. As written, the model 
has a Cu-0 hybridization €L = ±tpd with the minus sign occurring j = i+x/2 and j = i—y/2 
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and also hybridizaton = ±tpp between oxygen sites with the minus sign occurring for 
k = j— x/2~y/2 and k = j + x/2 + y/2. These phase conventions are illustrated in Fig. 1. 

The system is said to be half-filled when there is one hole per Cu site. At half-filling, for 
a wide range of parameters, its ground-state is an anti-ferromagnetic insulator just like the 
one-band Hubbard model. Its insulating gap , however, is a charge transfer one in contrast 
with the gap in the the one-band model which is a Mott-Hubbard gap. 

The values of the parameters in the Hamiltonian have been estimated by a number 
of different constrained density functional and quantum cluster calculations.0~0 Taken 
together, these estimates define a reasonably limited range of the parameters for the which 
the model might be labeled as "physical." . We use Table I of Ref 20, which summarizes 
results from density functional and cluster calculations, as the principal guideline for the 
parameter range we explored. 

The Cu site Coulomb repulsion Ud is large, making doubly occupancy of Cu sites by two 
holes very unfavorable. The next largest parameter, the charge transfer energy e = — > 
0, plays a special role. Depending on the relative values of e, Ud and the bandwidth W the 
system can be classified in different regimes .0 Estimates place the cuprate superconductors 
in the charge transfer regime, in which W < e < Ud- The remaining parameters are 
relatively comparable in magnitude. The role of the Cu-0 hybridization tpd is also special. 
This hybdization, through the super-exchange mechanism, generates an anti-ferromagnetic 
exchange interaction between the spins in the Cu sites. 

If e Ud, the three-band model maps into a one-band model with t^ff ~ ipd/e and 
U = Ud- For Ud ^ teff, the one-band model can in turn be mapped into the t-J model with 
J = Atljf/Ud- Zhang and Rice@ have argued that the t-J model can also be appropriate 
when tpd ^ e, Ud, Ud — e and e < Ud- In real materials, e/tpd is estimated to be ~ 2.7 — 3.7 .0 
Therefore, besides the lack of conclusive evidence that the one-band model superconducts, 
it is also unclear that the mapping among the most studied models is appropriate for phys- 
ical values of the parameters. In this paper we will examine the three-band model in the 
estimated range of parameters to detail the behavior of the model in its own right. 

5 



In what follows we scale all the energies by tpd- We take Vpd = since it is small (Vpd < 
O.StpJll). This assumption simplifies the numerical algorithm significantly and reduces the 
computational time by at least a factor of four. 

With the numerical method used, a variety of expectation values can be computed. We 
focused on two key quantities: the binding energy for doped holes and the pairing correlation 
functions as a function of the distance. The binding energy is defined as: 

A = (^2 - ^o) - 2 (El - Eo) (2) 

with En the ground-state energy for n doped holes, where Eq corresponds to the undoped 
(half-filled) case. If negative, this parameter indicates that it costs less to dope a second hole 
in the vicinity of the first hole than it does to create two isolated holes. In the thermodynamic 
limit A is expected to be the binding energy of a Copper pair. We calculated A for systems 
with 2x2,4x2 and 6x4 unit cells. We are unaware of any previous quantum Monte Carlo 
calculations of binding energy for the three-band model. 

We also computed the extended s-wave and the d^^-y^ pairing correlation functions: 

P,(i?) =< At (/?)A,(0) > (3) 

where 

<5 

with 6 = ±x, ±y. For extended s-wave pairing fs* (S) = 1 for all 6 and for the dx2_y2 pairing, 
/rf(5) = 1 for 5 = ±x and /d(5) = — 1 for 5 = ±y. Here R denotes the distance between 
unit cells (taken from copper ion to copper ion). We will directly report these correlations 
functions as a function of R. We will not report its g = spatial Fourier transformation 
and partial sums like Sa{L) = Y,r<l PaiR) as it has been done in previous works0^0, 
because the magnitude of these quantities are dominated by a large peak in Pa{R) when 
R is less than a few nearest neighbor distances. Over these distances. Pa measures local 



correlations among spin and charge fluctuations and has httle information about long-range 
pairing correlations. In certain cases, we will also report the "vertex contribution" to the 
correlation functions (see, for example. White et. aM) deflned as follows: 

Va{R) = Pa{R)- Pa {R) (5) 

where Pa (R) is the contribution of two dressed non-interacting propagators. For each term 
in Pa{R) of the form < cjc^cjcj > , Pa (R) has a term like < cjcj >< cjcj >. We flnd that 
the conclusions remain the same no matter which quantity we look at. 

III. CONSTRAINED PATH QUANTUM MONTE CARLO TECHNIQUE 

Our numerical method is extensively described and benchmarked elsewhere.@ Here we 
only discuss its basic approximation. In the CPMC method, the ground-state wave function 
|\l/o) is projected from a known initial wave function {"^t) by a branching random walk in 
an over-complete space of Slater determinants In such a space, we can write |\I^o) = 
J2<j,x{4')\4') ■ The random walk produces an ensemble of called random walkers, which 
represent |\l/o) in the sense that their distribution is a Monte Carlo sampling of that 
is, a sampling of the ground-state wave function. 

To completely specify the ground-state wave function, only determinants satisfying 
(^o|0) > are needed because |\l/o) resides in either of two degenerate halves of the Slater 
determinant space, separated by a nodal plane Af that is deflned by (\I/o|0) = 0. The sign 
problem occurs because walkers can cross A/" as their orbitals evolve continuously in the 
random walk. Asymptotically they populate the two halves equally, leading to an ensemble 
that has zero overlap with |\&o)- If A/" were known, we would simply constrain the random 
walk to one half of the space and obtain an exact solution of Schrodinger's equation. In 
the constrained-path QMC method, without a priori knowledge of A/", we use a trial wave 
function and require (\I/t|0) > 0. The random walk again solves Schrodinger's equation 
in determinant space, but under an approximate boundary-condition. This is what is called 
the constrained-path approximation. 
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The ground-state energy computed by the CPMC method is an upper bound. The quahty 
of the calculation clearly depends on the trial wave function {"^t)- Since the constraint only 
involves the overall sign of its overlap with any determinant it seems reasonable to 
expect the results to show some insensitivity to {"^t)- Through extensive benchmarking on 
the Hubbard model, it has been found that simple choices of this function can give very 
good results.! 

Besides as starting point and as a condition constraining a random walker, we also use 
|\E'r) as an importance function. Specifically we use (\E't|0) to bias the random walk into 
those parts of Slater determinant space that have a large overlap with the trial state. For 
all three uses of I^E^r), it clearly is advantageous to have {"^t) approximate |\l/o) as closely as 
possible. Only in the constraining of the path does |\E't) 7^ |\['o) generate an approximation. 

All the calculations reported here are done with periodic boundary conditions. Mostly, 
we study closed shell cases, for which the corresponding free-electron wave function is non- 
degenerate and is translationally invariant. In these cases, the free-electron wave function, 
represented by a single Slater determinant, is used as the trial wave function iV'r >• (The use 
of an unrestricted Hartree-Fock wave function as\ipT > produced no significant improvement 
in the results). 

However, to calculate the binding energy for holes, we need to study the half filled case 
and then the 1 and 2 hole doped cases. In the systems considered (2 x 2, 4 x 2 and 6x4 
unit cells) the 2 hole doped case corresponds to a closed shell case. The one hole doped case 
is a closed shell minus one hole. This missing hole can be put in either of two degenerate 
orbitals corresponding to ki = (0, vr) or k2 = (vr, 0). The trial wave function, constructed by 
putting the hole in either of these orbitals or a linear combination of them, give equally good 
results so we pick an arbitrary linear combination of them. For the half-filled case, the free 
electron wave function is 4-fold degenerate: ■^cl._^^^\CS >, c|.2_|c|,2 JCS* >, ^c^^ i\CS >, 
cl^^^cl_^^^\CS >, where \CS > is the closed shell state with two holes below half-filling. In 
the 2x2 system the state \CS > has 2 holes, in the 4x2 system, it has 6 holes and in 
the 6x4 system it has 22 holes. If we just arbitrarily pick one of these states, the energy 
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has a significant difference. For example, in the 2x2 system it is of the order of 10~^ from 
the exact diagonahzation result. To accurately compute the binding energy, this difference 
is too big. Therefore, we used the following procedure: we diagonalize the interacting 
part of the Hamiltonian in this degenerate subspace, and obtained 2 states with energy 
proportional to Ua and 2 states with zero energy. Of the 2 states with zero energy only 
one of them is a singlet. We used this state. It is represented by two Slater determinants: 
(^Ii,T^Ii,j, ~ ^L.T^L.J.) 1^*^ ^ / ■'■^ Table 1, for a 2 x 2 system with Cp — 3 and Up — tpp — 0, 
we compare the energies obtained using the CPMC with the one and two Slater determinants 
trial wave function and energies obtained from exact diagonahzation. We see that using two 
Slater determinants improves the accuracy by an order of magnitude or more. The accuracy 
has become better than the closed shell case. 

In a typical run for a large system we set the average number of walkers to 600 and the 
time step to 0.03. We perform 3000 steps before we start taking measurements and we do 
the measurements in 40 blocks of 500 steps each to ensure statistical independence. 

IV. RESULTS 

Shell structures are characteristic of finite-sized systems of electrons. In our simulations 
these effects are most clearly seen in the values of the energy and chemical potential as a 
function of the hole density n = N^/N, with N the number of unit cells. A shell structure is 
perhaps most easily illustrated by first considering the non-interacting problem. There, holes 
are added to the system at the same energy until all the degenerate states of a given shell are 
occupied and paired, that is, until the shell is closed. As an "open" shell is filled, the total 
energy of the system varies linearly with the doping and within a shell the chemical potential 
thus is a constant. The number of states in a shell is twice the order of the point group 
relating the wave vectors of the degenerate states. Between shells the chemical potential is 
discontinuous. This "gap" is a finite size effect. 

In the interacting problem the energy and chemical potential shows a similar behavior. 
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The number of states in a shell and the shell boundaries are the same as the non-interacting 
problem for the range of interactions we examined. For a 4 x 4 system, a typical shell 
structure for the three-band model as a function of Ud is illustrated in Figs. 2a and 2b were 
we show the ground state energy Eg and the chemical potential as a function of the hole 
density for e = 3 , tpp = t/p = 0. For our finite system we defined 

IJ.iN) = Eo{N)-EoiN-l) (6) 

which is the discrete version of the usual definition ^ = OEo/dNh- The Ud = ^ case 
is the non-interacting case. Its shell structure is most evident in the chemical potential, 
Fig. 2b. A comparison of the non- interacting case with the interacting cases does show 
several differences. At half-filling, n = 1, the chemical potential in the interacting cases 
shows a jump whereas the non-interacting case does not. This jump is not a finite-size 
effect but is the charge-transfer gap induced by the interaction. Within a shell the chemical 
potential of the interacting cases are varying weakly but approximately linearly with doping. 
This variation is similar to the behavior of the one-band Hubbard modeli'i, but there, to 
an excellent approximation, within a shell the energy for the interacting cases varies linearly 
with doping and hence within a shell the chemical potential is a constant. In Fig. 2a, 
one sees that the energy versus doping curve has a minimum whose location varies as a 
function of Ud- As Ud increases this location shifts from the fully doped position towards the 
half-filled position. The one-band Hubbard model displays similar doping and interaction 
dependencies. 

A. Binding Energy for Holes 

We calculated the binding energy A for holes for systems with 2 x 2, 4 x 2, and 6x4 
unit cells. A value of A < indicates hole binding whereas a value of A > implies hole 
repulsion. It is important to note that studying systems larger than 2 x 2 is essential for 
the doping 5 to have values similar to those for which the real materials superconduct: for 
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2x2 unit cells with one hole, S = 0.25, and with 2 holes, S = 0.5, values which are outside 
the physically interesting region. On the other hand, for one and two holes the 4x2 case 
has 5 = 0.125 and 5 = 0.25 while the 6 x 4 case has 5 = 0.042 and 5 = 0.083. These values 
of doping span the range observed in the real materials. 

Shown in Fig. 3, to benchmark the accuracy of our results, is a comparison of our binding 
energies with those of exact diagonalization. Specifically, we plotted A as a function of Ud 
for 2x2 systems with Up = tpp = and e = 2 and 3 and found good agreement with 
exact diagonalization results, especially for the e = 2 case. For e = 3 the differences are 
more pronounced but are still only of the order of 10^^. As a general trend, we slightly 
underestimate the binding energy. This figure is representative of our systematic error in 
the estimation of the binding energy. 

In Fig. 4 we present results for 4x2 systems as function of e as function of two values of 
Ud- For ?7rf = 1 we find that A decreases monotonically with e but for f/^ = 4 it has a shallow 
minimum at e ~ 2. Such a minimum means that there is a broad range of e for which the 
binding energy does not change significantly. The existence of such a minimum, occurring 
at approximately the same value of e, was observed in the exact diagonalization studies of 
Ogata and Shiba0 and commented upon by Martin.0 In general holes bind provided e is 
below some threshold value. In what follows, we set e = 3 which is a value consistent with 
the values obtained by LDA and cluster calculationjHSl and is also near the optimal value for 
binding. 

In Fig. 5, we compare A for different systems sizes as a function of Ud- We find that the 
binding energy tends to saturate as Ud is increased. The saturation is consistent with the 
exact diagonalization studies of Ogata and ShibJ and of Stephan et al.i The binding energy 
tends to increase as the system size is increased. Over the range of system sizes studied, 
the dependence of A on system size was not systematic enough to allow us to extrapolate 
to infinite system sizes. In general A ~ 10^^. 

We also studied the effect of Up and tpp on A. In Fig. 6 we plot A as a function of Ud 
for Up = and Up = 1- In Fig. 6a, we see that in the 2x2 system a small value of Up 
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destroys the binding as previously reported by Hirsch et al.Q However, in the larger 4x2 
system (Fig. 6b), the holes still bind with Up = 1, although the binding energy is smaller, 
confirming a speculation of Hirsch et al. that the negative effect of Up on the binding might 
decrease as the system size increases. A similar effect is seen when we plot A versus tpp for 
2x2 and 4x2 systems with f/^ = 6, e = 3 and Up = (Fig. 7). The ability to bind decreases 
nearly linearly with tpp but in the larger systems it decreases more slowly, so that there is 
still binding for tpp = 0.3. Therefore, the effect of the Coulomb repulsion on the O sites 
and hopping between them in reducing the binding is diminished for larger system sizes. In 
contrast to the position emphasized by Hirsch et al.i and by Stephan et al.,i the physical 
presence of a non-zero tpp and Up, while reducing the ability to bind, does not necessarily 
mean that an unphysically large Cu-0 Coulomb repulsion Vpd is essential for hole binding. 
Their conclusions were influenced by the small system sizes they were able to study. 

While the binding of two holes in the physically relevant range of parameters is encour- 
aging for superconductivity, perhaps a more central question now is what happens when 
a third hole is added. Do the three holes cluster or does the third hole sit apart from the 
pair? We do not examine these questions directly. Instead, we calculate the superconducting 
pairing correlations. 



B. Pairing Correlation Functions 

For both the extended s-wave and the c?^.2_y2 channels, we studied the pairing correlation 
functions as a function of the distance for lattices up to 6 x 6 unit cells. All calculations 
were done for closed shell cases. For the 4x4 systems, we put 18 holes {6 = 0.125), for 
the 6 X 6 system, 42 holes {6 = 0.167), and for 6 x 4, 26 holes {6 = 0.083). The distance is 
measured between unit cells (from copper ion to copper ion)ii. 

In Fig. 8 we present the pairing correlations for the extended s-wave channel of a 6 x 6 
system with e = 3 and tpp = 0.3 for different values of the interaction f/^. The inset magnifies 
the large distance behavior. We see that the short-ranged correlations are enhanced with 
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increasing Ud, but decay very quickly when R is increased. For large R they simply fluctuate 
around 0. For the same parameters, the d^2_y2 channel is plotted in Fig. 9. Again, there 
is enhancement in the short-ranged part. Then as R is increased, there is again a decay, in 
this case not to zero but to a small positive value. However, as Ud is increased, this value 
becomes even smaller. These results clarify the flndings of previous quantum Monte Carlo 
simulations. These previous simulations were limited to relatively high temperatures because 
of the sign problem. Scalettar and coworkers^ found an attractive pairing interaction in both 
the extended s and the dx'^-y'^ channels but could not predict which would dominate with 
the lowering of temperature and increasing of lattice size. The results of Figs. 8 and 9 clearly 
indicate that the dj.2_j^2 correlations dominate. This dominance is in direct opposition to the 
findings of Dopf and coworkers0'Ei that the extended s wave channel dominates. By looking 
at the the zero wave number component of the Fourier transform of the pairing correlation 
function, these researchers only observed how the short-ranged correlations behaved and 
overshadowed the behavior of the longer-ranged correlations. 

In several previous quantum Monte Carlo simulations various researchers found it more 
convenient to study just the vertex contribution to the pairing correlation function, as defined 
in section ||. In Fig. 10 we plot this contribution to the dx2_y2 correlation function just shown 
in Fig. 9. We see that there is positive contribution which again is enhanced by Ud at short 
distances but suppressed at larger distances. Thus, as emphasized by Zhang et al.,i this 
function seems to contain no more information about possible long-range order than the full 
pairing correlation function. 

Since we restricted ourselves to closed shell cases, it is difficult to do a useful finite- 
size scaling of our results. An attempt to show the effect of increasing system size on the 
correlations is given in Fig. 11 where we plot the dx2_y2 vertex contribution for 4x4 and 
6x6 systems with the same parameters, keeping in mind that the fillings are not the same. 
At short distances the correlations are enhanced as the lattice size is increased but for larger 
distances there is a crossover, and in fact, the vertex contribution is smaller in the 6x6 
system by almost an order of magnitude. A similar significant reduction is also seen in the 
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long-range parts of the full correlation functions. 

We also studied the effect of e and tpp on the pairing correlations. One will recall that 
increasing tpp weakened pair binding while the binding was optimal for e around 2. In Fig. 
12 we plot the bare correlation functions for different values of e with Ud — 6 and tpp — 0.3 
for a 6 X 4 system with 5 — 0.083. We find that the long-range pairing correlations are not 
optimal when the pair binding is optimal. In fact, values of e less than the optimal value 
enhance the longer range part of the correlation functions. A much smaller effect is obtained 
by changing the value of tpp, as evident in Fig. 13, where the parameters are Ud — 6, e = 3 
and the size is 6 x 6 unit cells. If anything, the larger value of tpp seems to enhance the 
long-range correlations slightly. 

It is interesting that in all cases presented we see a crossover from what happens at 
short distances to what happens at long distances. We emphasize that examining quantities 
that integrate, even partially, the pairing correlation function over distances can produce 
misleading results. 

V. SUMMARY AND CONCLUSIONS 

Using a newly developed quantum Monte Carlo method, the constrained-path method, 
we simulated the 3-band Hubbard model to search for a possible superconducting ground 
state. We focused principally on the calculation of the binding energy of a pair of holes 
and the extended s and dx2-.y2 superconducting pairing correlation functions. For the pair 
binding energy, we found that holes bind for a wide physically relevant range of parameters 
and that the binding increased as the system size is increased. We also found that when 
we included a hopping tpp between oxygen sites, the binding decreased roughly linearly with 
ipp 5 but in a larger system it decreased much more slowly. The Coulomb repulsion Up on 
the oxygen sites was found to suppress the binding, but again for larger systems its effect 
became less significant. The increased tendency to bind as the system size is increased seems 
to diminish the importance of an Cu-0 Coulomb repulsion Vpd for holes to bind as suggested 
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in several exact diagonalization studies-u'Ei In general, hole binding occurs over a broad range 
of e values centered roughly around two. In this respect our simulations support the exact 
diagonalization results of Ogata and Shiba.0 

We calculated the pairing correlation functions for the extended s and the dx2-y2 channels 
as a function of the distance. For the extended s channel, we found that as the Coulomb 
interaction Ud on the Cu sites is increased, the long-range part of the correlation functions is 
suppressed and fluctuates around zero. For the d^^-yi channel, we found that the correlations 
decay rapidly with distance towards a small positive value. However, this value becomes 
smaller as the interaction Ud or the system size is increased. The same behavior is observed 
for the vertex contributions to these correlation functions. Unequivocally, the dx2_y2 pairing 
correlations dominate the extended s pairing correlations. This finding is opposite to that 
previously emphasized by Dopf and coworkers.EHlil 

We also studied the effect of the charge transfer energy e on the pairing correlation 
functions and found that smaller values of e produce larger correlations. Similarly, we 
looked at the effect of the hopping tpp between oxygen sites and found that it has a very 
small effect but that larger values of tpp are slightly more favorable for pairing. These trends 
are counter to those that enhance pair binding. 

A systematic study of the size dependence of the binding energy and pairing correlation 
functions was difficult because the way the doping dependence of closed shell changed with 
changing system sizes. We did present a comparison of the d^2_y2 pairing correlation function 
for a 4 X 4 system and a 6 x 6 system. Although the hole doping was slightly different, a 
significant reduction of both the vertex contribution and the full correlation function was 
observed in the larger system. This result suggests that the long-range pairing correlation 
found are unlikely to persist in the thermodynamic limit. 

In general, we found similar trends for the three-band model that were found for the 
one-band Hubbard model.i The additional degrees of freedom in the three-band model do 
not seem to enhance superconductivity in an obvious way. One could argue that a Cu-0 
Coulomb repulsion Vpd is needed; however, we found that Vpd is not necessary to obtain hole 
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binding and speculate that it is unlikely that the addition of such a small parameter to the 
model would change the picture dramatically. The question remains, if the holes do bind, 
that is, a purely electronic mechanism induces an attractive interaction between holes, but 
this attraction does not lead to electron pairing, what does it, if anything, lead to. Does it 
lead to phase separation or some novel phase, like stripes? This question is currently under 
investigation. 
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TABLES 

TABLE I. Comparison of the exact ground-state energy with the CPMC result with a one 

and two Slater determinant trial wave function for a 2 x 2 system. Parameters are e = 3 and 
tpp = Up = 0. The use of two Slater determinants as described in the text improves the accuracy 
by an order of magnitude or more. 



Ud CPMC ISD CPMC 2SD Exact 

1 -5.0613(3) -5.0764(2) -5.076977 

2 -4.8475(9) -4.8789(7) -4.880047 
4 -4.6073(9) -4.6615(6) -4.661723 

6 -4.4884(9) -4.5468(6) -4.547436 
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FIGURES 

FIG. 1. Phase convention for the hopping matrix elements. The copper d^2_y2 orbital is 
surrounded by the oxygen and Py orbitals. The hopping matrix elements are shown with their 
corresponding phase. 

FIG. 2. (a) Ground state energy per unit cell as a function of the hole density for a 4 x 4 system 
with e = 3 , tpp = Up = 0. The energy decreases monotonically for Ud = but goes through a 
minimum for Ua > 0. For large Ud the minimum appears at half-filling, (b) Chemical potential as 
a function of the hole density for the same parameters. The shell structure is evident. For large 
Ud a charge transfer gap appears at half-filling. 

FIG. 3. Comparison of the binding energy as a function of Ud with exact diagonalization 

results for 2 X 2 systems with tpp = Up = and e = 2 and 3. Our results are in reasonable good 
agreement with the exact ones. 

FIG. 4. Binding energy A as a function of e for a 4 x 2 system with tpp = Up = 0. For Ud = 1 
it decreases monotonically but for [/d = 4 it has a shallow minimum near e ~ 2. 

FIG. 5. Binding energy A as a function of Ud for different system sizes with e = 3 and 
tpp = Up = 0. The binding increases as the system size is increased. 

FIG. 6. Effect of Up on the binding energy, e = 3 and tpp = 0. (a) In the 2x2 system the 
binding disappears completely for Up = 1. (b) In the 4x2 system the binding is still present for 
Up = 1 although is slightly suppressed. 

FIG. 7. Effect of tpp on the binding energy for 2 x 2 and 4x2 systems with e = 3, Ud = 6 and 
Up = 0. The binding energy decreases roughly linearly with tpp but it decreases much more slowly 
in the larger system. 
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FIG. 8. Extended s-wave pairing correlation functions vs. distance for different values of Ud 
for a 6 X 6 system with e = 3 and tpp = 0.3. The correlations decay rapidly and oscillate around 

zero at larger distances. 

FIG. 9. d-wave pairing correlation functions vs. distance for different values of Ud with the 
same parameters as in Fig. 8. The correlations decay towards a small positive value at large 
distances. 

FIG. 10. Vertex contribution for the d-wave pairing correlation functions vs. distance for 
the same parameters as in Fig. 8. It converges to a small value at large distances and this value 
decreases SiS Ud increases. 

FIG. 11. Comparison of the vertex contribution for the d-wave pairing correlations for two 
different system sizes with Ud = 6, e = 3 and tpp = 0.3. Although the fillings are slightly different 
in each case, a strong reduction is seen in the larger system. 

FIG. 12. d-wave pairing correlation functions for different values of e with Ud = 6,Up = tpp = 
for a 6 X 4 system with 6 = 0.083. The smaller values of e are more favorable for pairing. 

FIG. 13. Effect of tpp on the d-wave pairing correlation functions in a 6 x 6 system with Ud = 6 
and e = 3. A very small increase is seen as tpp is increased. 
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